Exploring Spectrum Sensitivies to Species Omission
For the sensitivity testing we estimated the size spectrum slope and intercept for every year and within each survey area an additional time, at each pass one of the species was omitted from the analysis. This was done to capture the relative importance each species had on the estimation of the biomass size spectra.
Loading Bio Data
The biological data used as a starting point is the same “node” used for the size spectrum estimations using the full community. This is loaded so that we can interrogate the overall biomass patterns that correspond to the results we see from the sensitivity testing.
---title: "Size Spectrum Species Sensitivity"author: name: "Adam A. Kemberling" title: "Quantitative Research Associate" url: https://github.com/adamkemberling affiliation: Gulf of Maine Research Institutedescription: | Supplemental materials for size spectrum manuscript.date: "`r Sys.Date()`"format: html: self-contained: true code-fold: true code-tools: true df-print: kable toc: true toc-depth: 3editor: sourceexecute: echo: false warning: false message: false fig.align: "center" comment: "" fig.width: 8 fig.height: 6---```{r setup}#### Packages ####library(tidyverse)library(gmRi)library(targets)library(ggridges)library(ggHoriPlot)library(patchwork)library(directlabels)library(treemapify)# set themetheme_set(theme_gmri())````r use_gmri_style_rmd(css_file = "gmri_rmarkdown.css")`# Exploring Spectrum Sensitivies to Species OmissionFor the sensitivity testing we estimated the size spectrum slope and intercept for every year and within each survey area an additional time, at each pass one of the species was omitted from the analysis. This was done to capture the relative importance each species had on the estimation of the biomass size spectra.## Loading Bio DataThe biological data used as a starting point is the same "node" used for the size spectrum estimations using the full community. This is loaded so that we can interrogate the overall biomass patterns that correspond to the results we see from the sensitivity testing.```{r}# Access Biological Data with {gmRi}# nmfs_bio <- gmRi::gmri_survdat_prep(survdat_source = "bio")withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(catch_1g_labelled))# Get Yearly and Regionally averaged Biomassesyr_szn_overall <- catch_1g_labelled %>%group_by(Year, survey_area) %>%summarise(overall_bio =sum(strat_total_lwbio_s, na.rm = T),overall_abund =sum(strat_total_abund_s, na.rm = T),.groups ="drop")# Get what it is for each speciessp_yr_szn <- catch_1g_labelled %>%group_by(comname, hare_group, fishery, Year, survey_area) %>%summarise(sp_bio =sum(strat_total_lwbio_s, na.rm = T),sp_abund =sum(strat_total_abund_s, na.rm = T),.groups ="drop")# Join in the overall, get the fraction and what the percent would beperc_bio <-left_join(yr_szn_overall, sp_yr_szn, by =c("Year", "survey_area")) %>%mutate(bio_frac = sp_bio / overall_bio,abund_frac = sp_abund / overall_abund,bio_perc = bio_frac *100,abund_perc = abund_frac *100)```### Overall Biomass Distributions```{r}# perc_bio overall rankingsregion_overalls <- perc_bio %>%group_by(comname, hare_group, fishery,survey_area) %>%summarise(avg_bio_perc =mean(bio_perc, na.rm = T),avg_abund_perc =mean(abund_perc, na.rm = T),.groups ="drop")```::: panel-tabset#### Biomass```{r}#| fig.height: 8# Treemapsggplot(region_overalls) +geom_treemap(aes(area = avg_bio_perc, fill = comname)) +geom_treemap_text(aes(area = avg_bio_perc, fill = comname, label = comname)) +facet_wrap(~survey_area, ncol =2) +scale_fill_gmri() +theme(legend.position ="none") +labs(title ="Fraction of Stratified Biomass: 1970-2019")```#### Abundance```{r}#| fig.height: 8# Treemapsggplot(region_overalls) +geom_treemap(aes(area = avg_abund_perc, fill = comname)) +geom_treemap_text(aes(area = avg_abund_perc, fill = comname, label = comname)) +facet_wrap(~survey_area, ncol =2) +scale_fill_gmri() +theme(legend.position ="none") +labs(title ="Fraction of Stratified Abundance: 1970-2019")```:::### Yearly Patterns::: panel-tabset#### Biomass```{r}#| fig.height: 8# Play around with plotting yearly influencesggplot(perc_bio, aes(Year, bio_perc)) +geom_col(aes(fill = hare_group)) +facet_wrap(~survey_area, ncol =1) +scale_color_gmri() +scale_fill_gmri() +theme(legend.position ="bottom") +labs(y ="Percent of Yearly Biomass",x =NULL,fill ="",title ="Percent of Stratified Biomass")# ```#### Abundance```{r}#| fig.height: 8# Play around with plotting yearly influencesggplot(perc_bio, aes(Year, abund_perc)) +geom_col(aes(fill = hare_group)) +facet_wrap(~survey_area, ncol =1) +#facet_grid(hare_group~survey_area) +scale_color_gmri() +scale_fill_gmri() +theme(legend.position ="bottom") +labs(y ="Percent of Yearly Abundance",x =NULL,fill ="",title ="Percent of Stratified Abundance")# ```:::## Loading Sensitivity Results```{r}# Access Biological Data with {gmRi}# nmfs_bio <- gmRi::gmri_survdat_prep(survdat_source = "bio")withr::with_dir(rprojroot::find_root('_targets.R'), tar_load(species_sensitivity_shifts))# Make year a numeric valuesens_results <- species_sensitivity_shifts %>%mutate(Year =as.numeric(Year)) # Join in the biomass info#names(sens_results)[which(names(sens_results) %in% names(perc_bio))]sens_results <- sens_results %>%select(-c(hare_group)) %>%left_join(perc_bio,by =c("Year"="Year", "survey_area"="survey_area", "spec_omit"="comname")) %>%drop_na() # All species are not present in all areas each year```### Single Species Removal Impacts```{r}# Get averages across all areas and years to rank themsens_ranks <- sens_results %>%group_by(`Species Removed`= spec_omit, Area = survey_area, Guild = hare_group, Fishery = fishery) %>%summarise(`Avg. Slope Shift`=mean(slope_shift, na.rm = T),`Avg. Slope Shift (z)`=mean(slope_shift_z, na.rm = T),.groups ="drop") %>%mutate(across(where(is.numeric), round, 2))```::: panel-tabset#### Slope Increase```{r}# Slope Steepenerssens_ranks %>%arrange(`Avg. Slope Shift (z)`) %>%head(8) %>% gt::gt() %>% gt::tab_header(title ="Top 8 Steepeners",subtitle ="Removing these species causes slope to steepen" )```#### Slope Flattener```{r}# Slope Flattenerssens_ranks %>%arrange(desc(`Avg. Slope Shift (z)`)) %>%head(8) %>% gt::gt() %>% gt::tab_header(title ="Top 8 Flatteners",subtitle ="Removing these species causes slope to flatten" )```:::### Regional Influence Heatmaps```{r}# Do we need to show the relative influence?# Things seem to wash out with the raw z-scores# Maybe they should be adjusted/weighted for how much any single omission# pulls a given year compared to either the # average z-shift for that year+region# or by the max that it shifts....``````{r}# plot heatmapinfluence_heatmap <-function(sens_results, filt_area, lims =c(-1, 1),rescale_z = F){# Make breaks to indicate limits are truncating blabs <-seq(from = lims[1], to = lims[2], length.out =5) blabs[c(1,5)] <-c(str_c(">", blabs[1]), str_c("<", blabs[5]))# Filter data filt_dat <- sens_results %>%mutate(Year =as.numeric(Year)) %>%filter(survey_area == filt_area)# Rescale?if(rescale_z == T){message("Rescaling not finished") }# Make heatmap filt_dat %>%ggplot(aes(Year, fct_rev(spec_omit), fill = slope_shift_z)) +geom_tile() +scale_fill_distiller(palette ="RdBu",limits = lims, labels = blabs,oob = scales::oob_squish, na.value ="gray20") +scale_x_continuous(expand =expansion(add =c(0,0))) +facet_wrap(~survey_area, ncol =2, scales ="free") +theme(legend.position ="bottom", legend.title =element_text(angle =0), axis.text.y =element_text(size =7),panel.background =element_rect(fill ="gray40"),panel.grid =element_line(color ="transparent")) +guides(fill =guide_colorbar(title.position ="top", title.hjust =0.5, barwidth =unit(3.25, "in"))) +labs(y ="Species omitted", fill ="Shift in Slope (z-score)")}```::: panel-tabset#### Gulf of Maine```{r}#| fig.height: 8# GOMsens_results %>%#filter(hare_group == "Groundfish") %>% influence_heatmap(filt_area ="GoM", lims =c(-2,2))```#### Georges Bank```{r}#| fig.height: 8# GOMinfluence_heatmap(sens_results, "GB", lims =c(-2,2))```#### Southern New England```{r}#| fig.height: 8# GOMinfluence_heatmap(sens_results, "SNE", lims =c(-4,4))```#### Mid-Atlantic Bight```{r}#| fig.height: 8# GOMinfluence_heatmap(sens_results, "MAB", lims =c(-4,4))```:::### Correlation to Biomass Fraction```{r}sens_results %>%ggplot(aes(bio_frac, slope_shift_z)) +geom_point(alpha =0.3) +geom_smooth(formula ="y ~ x", method ="lm") +labs(x ="Fraction of Biomass", y ="Shift in Spectra Slope (z)")```### Regional Influence Horizon PlotsI think this might help get around the "wash out" of the heatmap. The color scale should be longer, or potentially just not center on white...### Singe Species Checks: ```{r}# plot a species:# multiple panel plot function# Function for plotting a species that stands outplot_spectra_influence <-function(sens_results, filter_spec, filter_area){# Grab the relevant data filt_d <- sens_results %>%filter(spec_omit == filter_spec, survey_area == filter_area)# Make Plot Label p_label <-str_c(str_to_title(filter_spec), " Omission")# 1. Actual Value Shift# 1a. Slope splot <- filt_d %>%ggplot(aes(Year, l10_slope_strat)) +geom_line(aes(color ="All Species"), group =1) +geom_point(aes(color ="All Species")) +geom_line(aes(y = omit_slope_strat, color = p_label), group =1) +geom_point(aes(y = omit_slope_strat, color = p_label)) +facet_wrap(~survey_area, ncol =1) +scale_color_gmri() +theme(legend.position ="bottom") +labs(color =NULL,y ="Biomass Spectrum Slope",x =NULL)# 2. What the z-score shift is z_shift <- filt_d %>%ggplot(aes(Year, slope_shift_z)) +geom_line() +geom_point() +labs(y ="Slope Shift (z)",x =NULL)# Biomass Fraction Each Year: bfrac <- filt_d %>%ggplot(aes(Year, bio_frac)) +geom_line() +geom_point() +scale_y_continuous(labels = scales::percent) +labs(y ="Biomass Fraction")# Do they fit plotsemble <- splot / z_shift / bfracreturn(plotsemble)}```::: panel-tabset#### Haddock (GOM)```{r}#| fig.height: 8plot_spectra_influence(sens_results = sens_results, filter_spec ="haddock", filter_area ="GoM")```#### Spiny Dogfish (SNE)```{r}#| fig.height: 8plot_spectra_influence(sens_results = sens_results, filter_spec ="spiny dogfish", filter_area ="SNE")```#### Silver Hake (GoM)```{r}#| fig.height: 8plot_spectra_influence(sens_results = sens_results, filter_spec ="silver hake", filter_area ="GoM")```### Silver Hake (GB)```{r}#| fig.height: 8plot_spectra_influence(sens_results = sens_results, filter_spec ="silver hake", filter_area ="GB")```#### Butterfish (GOM)```{r}#| fig.height: 8plot_spectra_influence(sens_results = sens_results, filter_spec ="butterfish", filter_area ="GoM")```#### Atlantic Cod (GOM)```{r}#| fig.height: 8plot_spectra_influence(sens_results = sens_results, filter_spec ="atlantic cod", filter_area ="GoM")```:::
